COMMENT --------------------------------------------------------------------
This is a REDUCE program to calculate classical symmetries of ODEs

------------------------------- Initial data -------------------------------
 y(4) = a*y(2) + b*y(1) +c*y(0)$ is the initial ODE
 g is generating function of an arbitrary symmetry.
----------------------------------------------------------------------------;
OFF NAT$

OPERATOR x,y,g,a,b,c$
DEPEND a,x$ DEPEND b,x$ DEPEND c,x$

g:= v*y(1)+w$ 
DEPEND v,x,y(0)$ DEPEND w,x,y(0)$

%---------- Complete derivative operator of the initial ODE ----------------
FOR ALL ff LET 
CD(ff) = DF(ff,x) + y(1)*DF(ff,y(0)) + y(2)*DF(ff,y(1)) + y(3)*DF(ff,y(2))
       + (a*y(2) + b*y(1) +c*y(0))*DF(ff,y(3))$
%---------------------------------------------------------------------------
PROCEDURE CDK(f,k); % The k-power of CD.
BEGIN 
  SCALAR pp;
  pp:=f;
  IF k>0 THEN FOR i:=1:k DO pp:=CD(pp);
  RETURN pp;
END;
%---------- Universal linearisation operator of the initial ODE ------------
FOR ALL ff LET
LD(ff)=CDK(ff,4)-a*CDK(ff,2)-b*CD(ff)-c*ff$
%---------------------------------------------------------------------------
ORDER FACTOR y(3),y(2),y(1),y(0),x$

OUT result$
  g;
  LD(g); 
SHUT result$

END;